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Abstract Resource competition in heterogeneous en- 
vironments is still an unresolved problem of theoret- 
ical ecology. In this article I analyze competition be- 
tween two phytoplankton species in a deep water col- 
umn, where the distributions of main resources (light 
and a limiting nutrient) have opposing gradients and 
co-limitation by both resources causes a deep biomass 
maximum. Assuming that the species have a trade-off in 
resource requirements and the water column is weakly 
mixed, I apply the invasion threshold analysis (Ryabov 
and Blasius 2011) to determine relations between envi- 
ronmental conditions and phytoplankton composition. 
Although species deplete resources in the interior of 
the water column, the resource levels at the bottom 
and surface remain high. As a result, the slope of re- 
sources gradients becomes a new crucial factor which, 
rather than the local resource values, determines the 
outcome of competition. The value of resource gradi- 
ents nonlinearly depend on the density of consumers. 
This leads to complex relationships between environ- 
mental parameters and species composition. In partic- 
ular, it is shown that an increase of both the incident 
light intensity or bottom nutrient concentrations favors 
the best light competitors, while an increase of the tur- 
bulent mixing or background turbidity favors the best 
nutrient competitors. These results might be important 
for prediction of species composition in deep ocean. 
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1 Introduction 

Primary production forms the basis of metabolic ac- 
tivity of the ocean. Distinct phytoplankton groups con- 
tribute differently in the sequestration of CO2 (Frankig- 
noulle et al. 1994; Smetacek 1999), production of oxy- 
gen (Falkowski and Isozaki 2008), support of marine 
food webs (Christoffersen 1996), etc. A shift in the 
species composition may dramatically affect function- 
ing of the whole ecosystem (Walther et al. 2002; Cermeho 
et al. 2008; Paerl and Huisman 2009). However, in spite 
of the principal role of resource competition in the com- 
munity structuring, the conditions of coexistence and 
competitive exclusion in spatially variable environments 
still remain largely unknown. 

The classical theory, advanced by MacArthur (1972), 
Leon and Tumpson (1975), and Tilman (1980, 1982), 
analyses resource competition in uniform environments 
and shows that stationary coexistence of two species on 
two resources is possible only if growth of each species 
is finally restricted by its most limiting resource. The 
same results hold for competition in a mixed water col- 
umn where light exponentially decreases with depth 
(Huisman and Weissing 1995; Diehl 2002). However, 
in weakly-mixed systems these conditions may not be 
met. Competitors can be simultaneously limited by two 
or more resources, if their favorable habitats are sur- 
rounded by areas lacking these resources. 

For instance, in deep oligotrophic aquatic systems 
the light intensity reduces with depth, whereas concen- 
trations of nutrients typically have opposing gradients. 
As a consequence, species with distinct resource re- 
quirements can have maxima of production at different 
depths, which potentially decreases niche overlaps and 
increases biodiversity (Chesson 1990, 2000). However, a 
general extension of the competition theory to contin- 
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uous spatially variable habitats leads to difficult math- 
ematical problems (Grover 1997) and was addressed 
mainly in the mathematical literature (Hsu and Walt- 
man 1993). The analysis of this problem from the eco- 
logical point of view (Huisman et al. 1999; Yoshiyama 
et al. 2009; Dutkiewicz et al. 2009; Ryabov et al. 2010; 
Ryabov and Blasius 2011) is still far from complete and 
further research which would connect results for uni- 
form and spatially variable systems is required. 

In this article I analyze competition between two 
phytoplankton species in a deep weakly-mixed water 
column, assuming that limitation by light in deep lay- 
ers and limitation by nutrients at shallow depths cause 
deep chlorophyll or biomass maxima (Holm-Hansen and 
Hewes 2004; Kononen et al. 2003), which are a wide 
spread phenomenon in oligotrophic basins (Abbott et 
al. 1984; Karl and Letelier 2008). The location of a fa- 
vorable layer in such systems is not fixed, rather it de- 
pends on initial and boundary conditions, the stage of 
the relaxation process, etc. (Klausmeier and Litchman 
2001; Yoshiyama and Nakajima 2002; Ryabov et al. 
2010). Furthermore, a species, establishing at a certain 
depth, changes resource distributions and may affect all 
other species throughout the water column. Thereby 
this species acts as an ecosystem engineer, modifying 
its nutrient and light environment. 

For the analysis of competition in such a system 
Ryabov and Blasius (2011) recently introduced the no- 
tion of an invasion threshold, defined as a line (in case 
of two limiting resources) or a hypersurface (in gen- 
eral) in space of resource requirements, which separates 
the species that can grow in the presence of a resident 
species from those that cannot grow. The form and lo- 
cation of invasion thresholds depend on the characteris- 
tics of competitors as well as on the environmental con- 
ditions. The investigation of these dependences in the 
phytoplankton model reveals conditions that favor dif- 
ferent competitors and can explain shifts in the species 
composition caused by environmental changes. 



2 Model 

Competition between two phytoplankton species for light 
and a limiting nutrient (e.g., nitrogen or phosphorus) in 
a water column can be modeled in terms of a nonlocal 
system of coupled reaction-diffusion equations (Radach 
and Maier-Reimer 1975; Jamart et al. 1977; Klausmeier 
and Litchman 2001; Huisman et al. 2006) 
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where Pi(z,t) is the population density of the phyto- 
plankton species i at depth z and time t, (ii(I,N) is 
the growth rate, which depends on the local values of 
the light intensity, I(z,t), and nutrient concentration, 
N(z,t), rrii is the mortality rate, D is the turbulent 
diffusivity, and a 2 ; is the nutrient content of a phyto- 
plankton cell. 

The light intensity decreases with depth owing to 
the absorption of light by water and phytoplankton 
biomass (Kirk 1994) 
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where i,- n is the intensity of incident light, K^ g is the 
water turbidity and ki is the attenuation coefficient of 
phytoplankton cells. 

Assume that both resources are essential (von Liebig's 
law of minimum) and the resource limitation of growth 
can be parametrized by the Monod kinetics (Turpin 
1988), then the growth rate of species i follows 
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where /i ma x,i is the maximal growth rate, and -Hjv^ and 
Hj t i are the half-saturation constants, which define the 
species resources requirements. 

The phytoplankton cells cannot diffuse across the 
surface and bottom of the water column 
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The surface is also impenetrable for the nutrient, and 
the nutrient concentration at the bottom is constant 
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Fig.[TJshows typical equilibrium resource and biomass 
distributions in this model. The yellow shaded area in- 
dicates the production layer where both resources are 
available and the growth rate excesses mortality, n(N, I) > 
m. In a non- uniform system the total net growth on a 
favorable area should be large enough to compensate 
for losses into adjacent unfavorable areas (Ryabov and 
Blasius 2008). Denote the depths, which confine this 
area, as zjy if N limits species growth and zi if light 
limits species growth. The resource availability reaches 
at these depths the critical values 

— , I* = Hi t i , (5) 



N* = 



rrii 



at which growth equals mortality. The location of the 
production layer depends on the critical resource values 
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Fig. 1 (a) Distributions of the nutrient concentration (black), light intensity (gray), and phytoplankton biomass (green) in 
the model Eqn, IU41 (b) Distributions of the nutrient and light (solid lines) and exponential fitting (dashed lines) according to 
Eqn. \7\ and [5] plotted in the logarithmic scale. The yellow shaded area is the production layer, where the growth rate excesses 
mortality since N > N* and I > I*, (c) Zero net growth isoclines of species 1 and species 2 (green and red dashed lines, 
respectively); invasion threshold (blue lines) in the presence of species 1 as the resident under the conditions listed in Tabled 
The blue solid line was calculated numerically, by test of more than 5000 invaders with different half-saturation constants; the 
blue dashed line is a first order approximation with slope 71 = Ci,i/cn,i in the log-log scale CEa. 1151) 



and environmental parameters. In the present article, to 
focus on the influence of resource gradients, the param- 
eters are chosen to reproduce deep biomass maxima: 
the production layer is located in deep layers and the 
biomass density vanishes at the bottom and surface (see 
Table Q] for parameters). 

The invasion analysis is applied to determine out- 
comes of competition. Initially one species (the resi- 
dent) grows alone during a sufficiently long time and 
then another species (the invader) can grow from a 
small biomass density. Two species coexist if each of 
them can invade in the presence of its competitor. For 
the given model parameters the distributions of biomass 
and nutrients approach to equilibrium after approxi- 
mately 500 simulation days. To make sure that the res- 
ident is at equilibrium this time interval is increased up 
to 2000 days, after which an invader can grow and the 
system is simulated for further 18000 days to obtain the 
final competition outcome. 

The initial distribution of nutrients change linearly 
from at the surface to Nb at the bottom. The phyto- 
plankton species have initially uniform distribution of 
small density, P = 100 cells/m 3 . For the numerical so- 
lution the partial differential equations were discretized 
on a grid of 0.25 m. The resulting system of ordinary dif- 
ferential equations was solved by the CVODE package 
(http://www.netlib.org/ode) using the backward differ- 
entiation method. 



3 Competition in a spatially variable 
environment 

Assume that the competitors trade off in resource re- 
quirements: species 1 is a better nutrient competitor 
and species 2 is a better competitor for light (Fig. [lb) , 
and consider the invasion of species 2, assuming that 
species 1 has reached an equilibrium distribution. 

3.1 Invasion threshold 

Although the resident species depletes resources in the 
middle of the water column, the light intensity at the 
surface and nutrient concentrations at the bottom are 
still high (Fig. [Tbo, b). Due to a difference in resource 
requirements the maximum of production of species 2 
can be shifted towards, for instance, the nutrient sup- 
ply, and even a strong requirement for high nutrient 
concentrations can be compensated by adaptation to a 
low light intensity and vice versa. As a result, the inva- 
sion threshold takes the shape of a curve (the blue solid 
line in Fig. [lb is calculated by invasion of approximately 
5000 species with different half-saturation constants). 

The shape of invasion thresholds depends on the re- 
source distribution. However, to find this dependence 
we should find the principal eigenvalues corresponding 
to the boundary problem of Eqs.Q~r[2](Hsu and Waltman 
1993; Grover 1997; Ryabov and Blasius 2011), which 
can be solved in general only numerically (Troost et al. 
2005). However, a simple analysis can be performed for 
competitors with close resource requirements. Then the 
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maximum production of invaders occurs in the vicinity 
of the resident production layer (Fig. \Tjp) and it is pos- 
sible to express the growth rate of invaders in terms of 
the growth rate of the resident. 



3.2 Resource distribution 

Due to the diffusion of cells an essential part of the res- 
ident biomass is located outside the production layer. 
Consequently the resources below and above this layer 
are depleted to values which are even smaller than the 
critical values of the resident and I*. For instance, 
the nutrient concentration above the biomass maximum 
can be few orders of magnitude smaller than N* (FigJTJa). 
Therefore, the resident species shapes resource gradi- 
ents both within and outside of its production layer. 

The shape of resource distributions in this area will 
play a key role for the further analysis. It is commonly 
assumed, that the light intensity is exponentially dis- 
tributed, while nutrient concentrations can be fitted 
by a line (see the gray and black line, respectively, in 
Fig. [T^). However, the linear part of the nutrient pro- 
file is typically observed in the deep layers below the 
biomass maximum, where the light rather than the nu- 
trient limitation determines species growth. Similarly, 
above the production layer the growth is nutrient lim- 
ited and the net growth rate is negative independently 
of the shape of the nutrient distribution. The shape 
of resource distributions is crucial within the produc- 
tion layer, where the local net growth rate strongly 
correlates with local resource availability. As shown in 
Fig IB 3 ; m this area both the light and nutrient distribu- 
tions closely follow exponential dependences (straight 
dashed lines on the logarithmic scale). 

The exponential shape of resource distributions is 
a crucial point of the following theory. For this rea- 
son, it is important to note that this shape is not just 
an artifact of a specific model, but it was also found 
in field observations. For instance, Karl and Letelier 
(2008) clearly demonstrate exponential dependence of 
nutrient concentrations in the area of nutrient consump- 
tion. The emergence of such distributions can have dif- 
ferent nature. In particular, it is easy to show that it 
emerges, when the biomass variation within the pro- 
duction layer is small. 

Suppose that the phytoplankton density can be ap- 
proximated by a rectangular distribution with the con- 
stant density, 



Po = 



1 



P(z)dz 



within and in the vicinity of the production layer. This 
implies a small mortality level outside the production 



layer, m <C Umax, so that the biomass can diffuse with- 
out essential losses. Then, according to Eq. [31 the abso- 
lute value of the logarithmic gradient of light intensity 
is constant and equals 

d\nl(z) 



ci 



= K 



bg 



kP a . 



(6) 



Thus, the light distribution can be approximated as 

J{z) = I* e -cr(z-zt) . (7) 

To find the nutrient distribution note that accord- 
ing to Eqn. [5] in the limit of small mortality the criti- 
cal nutrient concentrations are also small, N* <C Hn, 
therefore the growth rate close to the depth zn can be 
linearized as jij^{N) » [i max N / Hn . Substituting this 
expression into Eq. [5] we obtain in equilibrium 
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A solution to this equation, that monotonically increases 
with depth, gives the equilibrium nutrient distribution 
in the vicinity of the depth zn, 

N{z)=N*e CN - l{z - ZN) 

with the logarithmic gradient 



(8) 
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CN = ^z— = y^i^r- (9) 

Thus in the vicinity of the production layer both the 
light and nutrient distribution have exponential shape, 
which implies that the logarithmic resource gradients 
have small variations with depth. 



3.3 Approximate calculation of the invasion thresholds 

To gain an insight into the dependence of invasion thresh- 
olds on resource distributions assume that the resident 
(species 1) and invader (species 2) differ only in their 
half-saturation constants, H^,i and Hj^, but are oth- 
erwise identical, i.e. fi max ,i = ^max,2, mi = m 2 , and 
D\ = -Da, (see Ryabov and Blasius (2011) for a general 
approach) . 

An invader of small initial density has a vanishing 
influence on the resources, therefore the possibility of 
invasion depends only on its growth rate in the resource 
distribution shaped by the resident, [i2(Ni(z), Ii(z)). 
Consider the difference between the nutrient limitations 
of the resident and invader. The Monod kinetics can be 
presented in the form 



N!(z)/H N . 
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which shows that the nutrient limitation of growth, fi^, 
depends only on the ratio N\/Hm,2- 
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Fig. 2 Schematic profiles of resource limitation of species growth, assuming that resource distributions are exponentially 
shaped, (a) If I(z) oc e~ c ' z (gray straight line on a logarithmic scale) then a difference in half-saturation constants of the 
species results in a parallel translation of the Monod function by Azi = —cj 1 In (Hi^/Hi^) along z-axis. In Figs.b-d solid 
lines show growth rates of each species, fi(N,I). The limitations of growth by particular resources, fipf(N) and fii(I), are 
shown as dashed and dot-dashed lines respectively in the regions where they do not coincide with fj,(N,I). The black dashed 
line shows the mortality rate, m. (b) Azi = Azn, parallel translation of the growth rate profiles, leads to the same total net 
growth of the species. In (c) and (d) either Azi < Azn or Azi > Azn and the net total growth rate of species 2 is either 
smaller or larger, respectively, than the total net growth of species 1. 
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However, the mathematical identity 

= , where Az = — In — 
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shows that division of an exponential function by differ- 
ent constants Hi and H^ is equivalent to a shift Az in 
position along the z-axis. Thus, using the exponential 
approximation for the nutrient distribution, Eq. [5J we 
obtain 



Ni(z) Ni(z-Az N ) 
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and the nutrient limitation of invader growth can be 
expressed through that of the resident 
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Therefore, the profiles of nutrient limitation of resident 
and invader growth have the same shape, but are shifted 
by Azn with respect to each other. If H^2 > Hn,i> 
i.e. the invader needs higher nutrient concentrations, 
then Azn is positive and the invader nutrient limita- 
tion profile is shifted downwards. Performing the same 
calculations for the light limitation we obtain 
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Azi has the opposite sign, because the light intensity 
has an inverse gradient. Assuming that the invader is 
a better light competitor (Hi 2 < Hm,i), we obtain 
that Azi is positive and the profile of light limitation 
is shifted downwards (see Fig. [5^). 

The values Azi and Azn show changes in the net 
growth arising from better adaptation to one or another 
resource, and the value A = Azn — Azi defines the dif- 
ference in the size of the resident and invader's favorable 
habitats. If boundary effects on biomass distributions 
are negligible, then the fate of the invader species de- 
pends solely on the sign of A. Indeed, if A = then 
the distinct resource requirements of these species re- 
sult in a parallel translation of the growth rate profile 
along the z-axis (Fig. [5b). Since the resident species has 
zero total growth in equilibrium, the same holds for the 
invader and its population cannot establish. If A > 
then the invader habitat is even smaller, consequently 
the invader total net growth will be negative (Fig. [5}i) . 
By contrast, if A < 0, the invader habitat is larger and 
it can invade the system because its total production is 
positive (Fig.[5H). 

According to Eqn. [5[ when the growth and mor- 
tality rates are equal, the ratio of half-saturation con- 
stants equals the ratio of critical resource values (e.g., 
Hi^/Hi i — I2/I1), and the invasibility threshold can 
be presented in the form 
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Dominance of sp 1 , y 



Fig. 3 (Left panel) Competition outcome in dependence of the slope of invasion thresholds of species 1 (green) and species 2 
(red). The gray dashed line shows the critical slope, ■y cr . (a) 71.2 < ~y C r - species 2 wins, (b) 72 < ~y C r < 7i - bistability, (c) 
71 < 7cr < 72 - coexistence, (d) 71,2 > 7cr - species 1 wins. (Right panel, e) All outcomes can be represented in dependence 
on 71 and 1/72- 



This inequality defines a first order approximation of 
the invasion threshold and has a straightforward geo- 
metrical interpretation: species 2 can invade in the pres- 
ence of species 1 if its critical resource values (iVj* , I|) 
lie below a line with slope 71 = ci_i/cn,i passing through 
the point (iV-*, /j") in the double- logarithmic resource 
plane (blue dashed line in Fig. [lb) and the range of pos- 
sible invaders is determined by the ratio of logarithmic 
resource gradients shaped by the resident. 

The value of 7 depends on the characteristics of 
the resident, and species 2 growing alone will shape a 
distinct resource distribution with different ratio 72 = 
ci.2/cn,2- Without loss of generality, assume that TV"! > 
TV-* and I| < A*i an d denote by 7 cr the slope (taken 
with opposite sign) of a straight line passing through 
the two critical resource points (gray dashed lines in 
Figs.[3ji-d), 

In/*//* 

Then, combining Eq. [TS] and its counterpart for species 
2, we obtain the outcomes of spatial resource competi- 
tion. If 71 < 7 CJ . < 72, then both species can invade 
the monoculture of each other and can thus coexist 
(Fig. [3b). In the opposite case, 71 > j cr > 72 none of 
the two species can invade, which leads to alternative 
stable states (Fig. [3b 1 ). Finally, one species dominates 
if 71,2 > lcr or 71,2 < 7cr (Fig. [3b, d). 

Consider in detail the conditions of coexistence. Fig. [3b 
shows that the best nutrient competitor (species 1, green) 
should have a relatively shallow slope (71 < j cr ) of its 
invasion threshold, therefore this species should shape 
a resource distribution with a relatively small gradient 
of light intensity, ci t \. This will provide more solar ra- 
diation for species 2 (red) , which, owing to its stronger 
nutrient limitation, has a niche in deeper layers. At the 



same time the invasion threshold of species 2 should 
have a steeper slope (72 > 7cr), and therefore cn,2 
should be relatively small. In other words, this species 
should not diminish the upward nutrient flux too much. 
Thus, we obtain a general rule that for coexistence each 
species should shape resource distributions with a rela- 
tively smaller gradient of its most limiting resource. 

It is convenient to represent all possible outcomes 
of competition in dependence of 71 and 1 /j2 (Fig. [3b) , 
which reflect the ability of a resident to shape a stronger 
gradient of its less limiting resource and quantify the 
dominance of the species over its competitor. Then co- 
existence is possible if the mutual dominance of both 
species is small, whereas large values of 71 and I/72 
lead to bistability. 

To get another perspective on the role of resource 
gradients, consider invasion in the presence of a resi- 
dent, which has shaped a large gradient of the nutri- 
ent distribution and a small gradient of light intensity, 
cn,i 3> ci t i. According to Eq. Q31 the change in the 
area of light limitation, Azi, approaches infinity when 
c/,i — > 0, thus even a small difference in light require- 
ment, ln(/|//*) < 0, can lead to a large increase of the 
favorable area. Therefore, a better adaptation to this 
resource is very efficient. By contrast, a large nutrient 
gradient cm,i makes competition much harder, because 
the areas of the resident and invader nutrient limita- 
tion will almost coincide (Azn oc 1/cn.i — > 0, unless 

N£ < N * or N 2 > N i, see E q- US- In thc limit case 
cn,i — > 00, competition for the nutrient becomes impos- 
sible because the habitat is virtually divided into two 
parts: with very low nutrient concentrations (all species 
are nutrient limited) and very high concentrations (the 
nutrient limitation plays no role). 



Phytoplankton competition in deep biomass maximum 



7 



N = 0.5 mmnl/m 

2 1 

Light, umol photons m s 

0,1 1 10 100 1000 
T" 



N = 1 mmol/m 

2 1 

Light, umol photons m s 

0,1 1 10 100 1000 
T 



N = 10 mmol/m 

2 1 

Light, umol photons m s 

0,1 1 10 100 1000 

r 




10 10 10 

Nutrients, mmol/m 



0,0 0,2 0,4 0,6 

Plankton, 10 7 cells m" ; 



10 10 10 

Nutrients, mmol/m 3 



0,0 0,2 0,4 0,6 0,8 1,0 1,2 

Plankton, 10 7 cells m~ 3 



N* N* 



Iff* 10" ! 10" 1 10° 10' 

Nutrients, mmol/m 3 

6 5 10 15 20 25 
Plankton, 10 7 cells m" 3 




Fig. 4 Geometrical method to project the outcome of spatial resource competition in dependence of the nutrient concen- 
tration at the bottom, Nb- Top panel: shift in the competition outcome between species 1 (green) and species 2 (red) in the 
phytoplankton model caused by an increase of Nb- Bottom panel illustrates this shift as the result of a counter-clockwise 
rotation of the invasion thresholds. The slopes, 71 and 72, of the invasion thresholds in the bottom panel were calculated 
numerically for a monoculture of species 1 and 2. 



4 The influence of environmental parameters 
on the competition outcome 

In this section the invasion threshold analysis is applied 
to explain shifts in the species composition predicted by 
model Eqn. HHU First two examples demonstrate how 
changes in the competition outcome can be visualized 
as rotation of invasion thresholds and then changes in 
the species composition are considered in the (NB,Iin) 
and (D,Kb g ) planes. 

4.1 Shifts in the species composition with D and Nb 

Compare changes in the composition of phytoplank- 
ton species, caused by an increase of turbulent mixing, 
D, and bottom nutrient concentration, Nb- Both these 
parameters increase the nutrient availability through- 
out the water column. In a homogeneous environment 
higher nutrient concentrations would provide better con- 
ditions for the better light competitor (species 2). How- 
ever, surprisingly, in the phytoplankton model these 
changes have opposite effect. While an increase of Nb 
favors the better light competitor (Fig. HJ), an increase 



of D favors to the better nutrient competitor (Fig. [SJ. 
To explain these differences consider the influence of 
these parameters on the resource gradients and the slope 
of invasion thresholds. 

More nutrients at the bottom give rise to a larger 
and sharper biomass distribution of the resident species, 
which in turn yields a steeper gradient of the nutrient 
concentration. As a result, with an increase of Nb, the 
low light adapted species 2 "shades" the nutrient flux 
more strongly and hinders invasion of a better nutri- 
ent competitor species 1, which, in consequence of its 
light limitation, occupies higher layers. Because a larger 
nutrient gradient cn.i leads to smaller values of 7*, this 
transition can be represented as a counter-clockwise ro- 
tation of the invasion thresholds in the resource plane 
(the bottom panel in Fig. 0]). 

By contrast, with an increase of mixing, D, the 
biomass distribution becomes wider and the nutrient 
gradient decreases (Fig. [5]). As a result, more nutri- 
ents become available for species 1, which finally wins 
the competition. In the resource plane this change can 
be represented as a clockwise rotation of the invasion 
thresholds. 
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D = 0.1 cm 2 /s D = 0.9 cm 2 /s D = 2.5 cm 2 /s 

Light, ^mol photons Light, ^mol photons Light, |j.mol photons 




AT* '° N* '° N' 

Fig. 5 The same as in Fig. ^but in relation to diflusivity. JVb = 2 mmol nutrient m -3 



4.2 Numerical simulations 

The value 7 quantifies the likelihood of invasion and can 
be used derive the competition outcome from the results 
obtained for each competitor alone. Consider species 1 
(the better nutrient competitor) as the resident. Fig. [6^, 
shows the values of arctan7i in the (Ji„,7Vs) plane, 
which were numerically calculated from equilibrium re- 
source distributions in the presence of this species. The 
intensity of red color decreases with 71 and character- 
izes the ability of a better light competitor (species 2) 
to invade. Contrary, the dominance of species 1 over 
species 2 increases with 7. The black line shows the 
level 71 = 7 cr . In the region to the right of this line 
7i < 7cr and therefore, species 2 can invade. In a sim- 
ilar manner, Fig. ^jp shows the slope of the invasion 
threshold, arctan72, in the presence of species 2 alone. 
Here, however, the color intensity increases with 72 and 
characterizes more favorable conditions for a better nu- 
trient competitor. The black line is the boundary of 
the range 72 > 7 cr , where species 1 can invade. An 
intersection of these ranges shows the range of coexis- 
tence, in which 71 < 7 cr < 72 and each species can in- 
vade in the presence of its competitor. Figj6j; compares 
this range, based on the one-species modeling (between 
black lines) , with the range of coexistence obtained by 
two-species modeling (the blue area). 



This approach can be extended to include other 
model parameters. For instance, Figs. [SJi-f show that 
an increase of both the turbulent diffusivity, D, and 
background turbidity, K\, g increases the slopes of the 
invasion thresholds and the competition outcome shifts 
from the dominance of species 2 through the range of 
coexistence to the dominance of species 1. In this figure 
again, the range of species coexistence predicted from 
simulation of one-species model (the area between two 
black lines in Fig. [Bf ) is in a good agreement with the 
results obtained in two-species model (the blue area). 

5 Discussion 

5.1 Invasion thresholds 

In a uniform system a species can increase its biomass 
if its resource requirements are lower than the present 
level of ambient resources (i?*-rule). This rule, how- 
ever, has to be generalized for spatially variable envi- 
ronments, where on the one hand the size of favorable 
area becomes a crucial factor, and on the other hand 
the resource heterogeneity provides an opportunity to 
compensate a lack of one resource by superfluous con- 
centrations of another. To extend the competition the- 
ory to nonuniform systems Ryabov and Blasius (2011) 
recently suggested to replace the i?*-rule by the notion 
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Fig. 6 Invasion analysis in the phytoplankton model. (Top panel) The slope of the invasion thresholds as a function of 
Nb and Ii n numerically calculated for a monoculture of (a) species 1 and (b) species 2. The color intensity changes with 
arctan ji and characterizes the likelihood of invasion of (a) better light competitors (decreases with 7) and (b) better nutrient 
competitors (increases with 7). Black lines shows the boundary of the parameter range where the competitor can invade, (c) 
Compares the boundary of the coexistence range obtained from results for monocultures (Figs, a and b) and from two-species 
modeling. The color shows the logarithm of the biomass ratio (logBi/B2) after 20000 simulation days. (Bottom panel) The 
same analysis in dependence on D and K^ g . 



of an invasion threshold, which is defined as the maxi- 
mal requirements of successful invaders in the presences 
of a resident species in equilibrium. If the critical val- 
ues (N*,I*) of a species lie below the threshold, then 
the species can invade. In a nonuniform environment 
the shape of invasion thresholds can be complex and 
nonlinear. However, an approximated technique can be 
developed for competition between species with close 
resource requirements. For these species the invasion 
threshold can be approximated by a strait line on a 
double-logarithmic scale and the slope of this line (de- 
termined by the ratio of logarithmic resource gradients) 
becomes a critical determinant of the competition out- 
come. 

The dependence of invasion thresholds on resource 
gradients, rather than on the local resource values, leads 
to new rules for invasion and coexistence. In particular 
a large value of 7 means that the invasion threshold 
approaches to a vertical line, which favors to good nu- 
trient competitors (see Fig. [3]). By contrast, if 7 is small 
then the invasion thresholds is close to a horizontal line, 
therefore good light competitors are in more favorable 
conditions. 



This effect alters the mechanism of coexistence. In a 
uniform system two species can coexist if each of them 
mostly consumes its most limiting resource, and finally 
becomes self-limited by this resource. If however, the 
favorable area is bounded by resource availabilities, a 
somewhat opposite rule can be formulated: the mono- 
culture of each species should shape resource distribu- 
tions with a relatively smaller gradient of its most limit- 
ing resource. Then each competitor may benefit from its 
adaptation to a specific resource. As shown in (Ryabov 
and Blasius 2011), for phytoplankton competition this 
difference is even more striking, because 7 grows with 
light attenuation coefficient, fc, and decrease with nutri- 
ent content of cells, a. Thus, for coexistence in deep lay- 
ers of a water column a low- light/high- nutrient adapted 
species should have a smaller value of a, so that the 
nutrient is available also for species in the upper lay- 
ers. Conversely, the species with high- light/low- nutrient 
adaptation should have a smaller coefficient k, thereby 
minimizing the light shading. Note that similar corre- 
lations between consumption rates and resource adap- 
tation were recently found in the experimental analysis 
of competition for light and phosphorus (Passarge et al. 
2006). 
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The dependence of competition outcome on envi- 
ronmental parameters in the phytoplankton model also 
reveals a number of intriguing results and shows that 
the intuition based on homogeneous models may fail in 
analysis of heterogeneous systems. For instance, Figs. [6k.- 
c show that an increase of both the incident light in- 
tensity, Ii n , and nutrient concentrations at the bottom, 
Nb, decreases 7 and favors to the best light competi- 
tor (species 2). This is because both factors lead to 
a sharper biomass maximum and hence to a steeper 
nutrient gradient. Recall that in a uniform system an 
increase of / would always favor to a good nutrient com- 
petitor, while an increase of N would favor to a good 
light competitor. 

Further, Figs. [SJi-f show that an increase of the wa- 
ter turbidity, Kb g , or turbulent diffusivity, D, increases 
7 and favors to the best nutrient competitor (species 1). 
Although leading to the same results, these changes in- 
volve different mechanisms: an increase of Kb g increases 
the ratio 7 = cj /cm via an increase of the light gradient, 
c/, whereas an increase of D leads to the same result 
because the nutrient gradient cjv becomes smaller. 

These results might have important outcomes for 
field research. For instance, a decrease of Nb or D de- 
creases the nutrient availability in the water column, 
however, leads to opposite effects on the competition 
outcome. Thus, a shift in the species composition caused 
by higher stratification of the ocean waters can be op- 
posite to that caused by a reduction of nutrient levels in 
deep layers. This effect can possibly explain both pos- 
itive and negative correlations between nutrient con- 
centrations and the abundance of high-light adapted 
species observed along environmental gradients in the 
Atlantic Ocean (Johnson et al. 2006). 

5.2 Model assumptions 

To derive the invasibility criterion (Eg. IT5|) . an "ideal" 
system was considered. It was assumed that the re- 
sources are exponentially distributed, the biomass max- 
imum is located in the deep layer and /i roax ,i = Mmax,2- 
Under these assumptions the analysis of competition 
becomes fairly simple, because the invasion threshold 
follows a straight line with slope 7 in double logarith- 
mic space. Now consider a more general situation in 
which these assumptions do not hold. 

First, if the resource distributions are not expo- 
nentially shaped then the invasion threshold takes the 
shape of a curve (blue solid line in Fig. QJ). However, 
a tangent line to this curve at the point (N*, Ij") is 
exactly the linearly approximated invasion threshold 
(blue dashed line). Thus, if resources deviate from ex- 
ponential distributions, Eq. \T5\ provides a first order 



approximation which is valid for species with close re- 
source requirements. The larger the interval where re- 
sources change exponentially, the larger the segment of 
the invasion threshold which follows the linear depen- 
dence in log-log scale. Often the resource level changes 
exponentially for few and more orders of magnitude 
(Ryabov & Blasius 2011). In this case the invasibility 
criterion Eq. [15] is applicable, if differences in species 
half-saturation constants have the same or smaller or- 
der. 

Second, it was assumed that the production layers 
are confined by the resource availability and located 
sufficiently far from the boundaries. The boundaries ef- 
fect species distributions and survival conditions. For 
instance, an impenetrable boundary is a favorable fac- 
tor for species survival, because the cell diffusive flux re- 
flects from the boundary and more cells can return into 
the production layer (Cantrell and Cosner 1991). The 
boundaries also influence the species spatial segrega- 
tion. The species, which can occupy different depths in 
the deep layers, have to compete locally if the biomass 
maximum occurs in the benthic or surface layer, which 
strengthens the interspecific interactions. All these ef- 
fects have a profound impact on species composition 
and can even reverse the outcome of competition, will 
be published elsewhere. 

Third, to focus on the role of resource gradients 
it was assumed that competitors have the same maxi- 
mal growth rates, mortalities and dispersal. In this set- 
tings an invasion threshold always passes through the 
resident's (N* , I*) values. A more general approach 
(Ryabov and Blasius 2011) shows, however, that if there 
is any difference in these parameters then the invasion 
threshold can be shifted towards higher or lower re- 
source requirements. As a result an invader with, for 
instance, higher fi max can invade even if it has higher re- 
quirements of both resources and two species can stably 
coexist via a positive correlation in the maximal growth 
rate and resource requirements (so-called gleaner-opportunist 
trade-off). 



5.3 Comparison with other models 

There are few approaches suggested to describe phy- 
toplankton competition in a water column. These ap- 
proached can be classified based on the assumed in- 
tensity of water mixing. Namely, Huisman and Weiss- 
ing (1994, 1995) consider competition in well-mixed 
systems; by contrast, Yoshiyama et al. (2009) suggest 
an approach for poorly mixed environments; finally, 
Ryabov and Blasius (2011) and the present paper com- 
plement these studies for the case of moderate mixing. 
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Compare the main assumptions and outcomes of these 
models. 

Huisman and Weissing (1994, 1995) performed an 
analysis, assuming that the light intensity decays ex- 
ponentially, but phytoplankton biomass and nutrients 
are uniformly distributed (Fig. [7^). The distributions 
of the competitors in this system completely overlap. 
The outcome of their competition changes with the ra- 
tio and absolute values of resource supplies. If com- 
petitors trade off in resource requirements, then (sim- 
ilar to Tilman 1980) the outcome of competition de- 
pends on the resource ratio. If however, one competitor 
has a higher /z max then this species can benefit both 
from light and nutrients. Moreover, even having higher 
requirements for both resources, this species can win 
competition if the resource supplies are high enough. 

Yoshiyama et al. (2009) consider competition in a 
stratified water column with a uniformly mixed upper 
layer and poorly mixed deep layer. For the upper layer 
the approach of Huisman and Weissing is applied, for 
the deep layer it is assumed that the biomass maximum 
is infinitely thin (Fig. [T]b) . According to the last as- 
sumption, if species compete in the deep layer their dis- 
tributions will never overlap. Thus, this model is more 
suitable if traits of two species are sufficiently different. 
For such species the model predicts that in the deep 
layer the resource supplies do not change the outcome 
of competition, they rather influence the depth of the 
biomass maxima. The outcome of competition changes, 
however, if the bulk biomass of at least one of the com- 
petitors occurs in the upper or bcnthic layer. 

The approach of Ryabov and Blasius (2011) consid- 
ers competition in the deep layer of a moderately or 
poorly mixed water column (Fig. [7£). This approach is 
based on the comparison of the growth rate profiles in 
the presence of each competitor alone (Fig. ^ . Conse- 
quently, the analysis does not depend on the biomass 
distribution, and the biomass maxima can overlap. Fur- 
thermore, this approach reveals a key role of resource 
gradients in community structuring. 

A single mathematical model cannot present exactly 
the dynamics of a complex ecological system. All mod- 
els contain some simplifications, and typically a real 
system and model match qualitatively, but not quanti- 
tatively. In this sense, the invasion thresholds can be- 
come a useful tool of the qualitative analysis. In the 
model considered here the invasion thresholds have a 
simple linear shape in a log-log scale. In another model 
this shape can be non-linear or linear, but in another 
scale. If, however, this shape and its dependence on 
system parameters can be deduced numerically or ana- 
lytically, than we can also project shifts in the species 
composition for this system. Thus, this analysis can be 



further extended for a wide spectrum of spatially het- 
erogeneous models, in which other biotic or abiotic fac- 
tors, such as gradients of temperature, predation, mor- 
tality, etc., are taken into account. 

6 Acknowledgments 

I thank B. Blasius for comments and useful suggestions. 



References 

1. Abbott MR, Denman KL, Powell TM, Richerson PJ, 
Richards RC, Goldman CR (1984) Mixing and the Dy- 
namics of the Deep Chlorophyll Maximum in Lake Tahoe. 
Limnol. Oceanogr. 29:862-878 

2. Cermeno P, Dutkiewicz S, Harris RP, Follows M, Schofield 
O, Falkowski PG (2008) The role of nutricline depth 
in regulating the ocean carbon cycle. Proceedings of 
the National Academy of Sciences 105:20344-20349. DOI: 
10.1073/pnas.0811302106 

3. Chesson P (1990) MacArthur's consumer-resource model. 
Theor. Popul. Biol. 37:26-38 

4. Chesson P (2000) Mechanisms of Maintenance of Species 
Diversity. Annual Review of Ecology and Systematics 
31:343-366 

5. Cantrell RS, Cosner C (1991). The effects of spatial hetero- 
geneity in population dynamics. J Math Biol, 29:315-338. 

6. Christoffersen K (1996) Ecological implications of 
cyanobacterial toxins in aquatic food webs. Phycologia 
35:42-50 

7. Diehl S (2002) Phytoplankton, Light, and Nutrients in a 
Gradient of Mixing Depths: Theory. Ecology 83:386-398 

8. Dutkiewicz S, Follows MJ, Bragg JG (2009) Modeling the 
coupling of ocean ecology and biogeochemistry. Global 
Biogeochem. Cycles 23:GB4017 

9. Falkowski PG, Isozaki Y (2008) The story of 02. Science 
322:540-542 

10. Frankignoulle M, Canon C, Gattuso JP (1994) Ma- 
rine calcification as a source of carbon dioxide: Positive 
feedback of increasing atmospheric CO2. Limnology and 
Oceanography :45 8-462 

11. Grover JP (1997) Resource Competition Springer 

12. Holm-Hansen O, Hewes CD (2004) Deep chlorophyll- 
a maxima (DCMs) in Antarctic waters. Polar Biology 
27:699-710. DOI: 10.1007/s00300-004-0641-l 

13. Hsu SB, Waltman P (1993) On a system of reaction- 
diffusion equations arising from competition in an un- 
stirred chemostat. SIAM J. Appl. Math 53:1026-1044 

14. Huisman J, van Oostveen P, F. J. Weissing (1999) Species 
Dynamics in Phytoplankton Blooms: Incomplete Mixing 
and Competition for Light. Am. Nat. 154:46-68 

15. Huisman J, Weissing FJ (1994) Light-Limited Growth 
and Competition for Light in Well-Mixed Aquatic Envi- 
ronments: An Elementary Model. Ecology. 75:507-520 

16. Huisman J, Weissing FJ (1995) Competition for Nutri- 
ents and Light in a Mixed Water Column: A Theoretical 
Analysis. Am. Nat. 146:536-564 

17. Huisman J, Pham Thi NN, Karl DM, Sommeijer B (2006) 
Reduced mixing generates oscillations and chaos in the 
oceanic deep chlorophyll maximum. Nature 439:322-325. 
DOI: 10.1038/nature04245 



12 



Alexei B. Ryabov 



(a) Unifomly mixed water column (b1) Upper mixed layer (UML) (b2) Poorly mixed deep layer (c) Moderately mixed deep layer 



Nutrients 



Nutrients 



Nutrients 



Nutrients 



N* 



TV* 



1 1 




1 1 1 


— i 
5 








Production 










layer 




















i i i 







Production 
layer 






• 





/* 

Light 





Light 



N* 







\ Production* — ^ 
\ layer^ — ' 


[ / 
i 




/* 

Light 



p 

-N 
I 



Fig. 7 Different modeling approaches for competition in a water column, (a) Huisman and Weissing (1994, 1995) assume an 
exponential decay of light and uniform distribution of nutrients and phytoplankton. The light limitation bounds the production 
layer (yellow area). Yoshiyama et al. (2009) assume (bl) within the UML the biomass and nutrient are uniformly distributed 
and (b2) in the poorly mixed deep layer the biomass distribution is infinitely thin. The production layer is bounded by the 
availability of light in the UML and by both resources in the deep layer, (c) Ryabov and Blasius (2011) assume exponential 
resource distributions. The production layer is bounded by the availability of nutrients (upper boundary) and light (lower 
boundary) and has finite thickness. The biomass distribution can have an arbitrary shape. 
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Table 1 Parameters values and their meaning 



Symbol 


Interpretation 


Units 


Value 




Independent variables 








t 


Time 


h 


- 




z 


Depth 


m 


- 




Dependent variables 








P(z,t) 


Population density 


cells m -3 






I(z,t) 


Light intensity 


/imol photons m -2 s -1 






N(z,t) 


Nutrient concentration 


mmol nutrient m -3 






Parameters 








/in 


Incident light intensity 


/imol photons m -2 s _1 


600 (100 - 


1000) 




Background turbidity 


m" 1 


0.1 




Z B 


Depth of the water column 


m 


100 




N b 


Nutrient concentration at Zb 


mmol nutrient m — 3 


2 (0.1-10) 




D 


Vertical turbulent diffusivity 


cm 2 s — 1 


0.3 (0.02 - 


5) 




Maximum specific growth rate 


h" 1 


0.04 




m 


Specific loss rate 


h" 1 


0.01 




Species 


1 — best nutrient competitor 








Hi.i 


Half-saturation constant, light limitation 


/imol photons m -2 s _1 


20 




Hn,i 


Half-saturation constant, nutrient limitation 


mmol nutrient m -3 


0.04 




ai 


Nutrient content of phytoplankton 


mmol nutrient cell -1 


8 xlO -10 




ki 


Absorption coefficient of a phytoplankton cell 


m 2 cell -1 


6xl0 -10 




Species 


2 — best light competitor 








Hj.2 


Half-saturation, light limitation 


/imol photons m - 2 s -1 


15 




Hn,2 


Half-saturation constant nutrient limitation 


mmol nutrient m -3 


0.065 




Ct2 


Nutrient content of phytoplankton 


mmol nutrient cell -1 


5 xlO -10 




k 2 


Absorption coefficient of a phytoplankton cell 


m 2 cell -1 


6xl0 -10 





